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A new quark-field smearing algorithm is defined which enables efficient calculations of a broad 
range of hadron correlation functions. The technique applies a low-rank operator to define smooth 
fields that are to be used in hadron creation operators. The resulting space of smooth fields is small 
enough that all elements of the reduced quark propagator can be computed exactly at reasonable 
computational cost. Correlations between arbitrary sources, including multi- hadron operators can 
be computed a posteriori without requiring new lattice Dirac operator inversions. The method is 
tested on realistic lattice sizes with light dynamical quarks. 
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I. INTRODUCTION 

One of the goals of lattice calculations is to predict the low-energy hadron spectrum of confined quarks and gluons, 
starting solely from the QCD lagrangian. This approach to spectroscopy necessitates methods for measuring the 
two-point correlation functions of field operators with the selected quantum numbers under investigation. A complete 
understanding of the QCD spectrum must also included excitations of the mesons and baryons as well as exotic 
states. Describing the resonances seen in scattering experiments as states in QCD has long presented a challenge to 
lattice calculations as direct access to the matrix elements related to decay widths is usually missing in Euclidean 
formulations of field theory. The way to circumvent this is well known in principle: decay properties can be inferred 
from a detailed study of the dependence of the spectrum in a finite box on the volume of the box [l|, Q ■ The best means 
of making solid determinations of energy levels is to employ a large basis of operators and then use the variational 
method to find good approximations to energy cigenstatcs. As soon as states above threshold arc to be explored, the 
basis should include operators that resemble multi-hadron systems where the constituents have well-defined momenta. 
Access to all elements of the quark propagator 0] relevant to long-range correlation enables these measurements and 
also enables more detailed investigations into physical observables whose accurate measurement has often eluded 
lattice practitioners. One example is the isoscalar meson sector, where mass determinations require the evaluation of 
disconnected diagrams as part of the two-point correlation function. As well as this ambitious list of requirements, 
precision data will be crucial for these calculations 0]. 

In a Monte Carlo calculation on the lattice, the physically relevant signal in a correlation function falls exponentially 
and is rapidly dominated by statistical fluctuations. Operators that create low- lying energy eigenstates quickly are 
invaluable and improve the quality of data extracted exponentially. The most useful tool at the lattice practitioner's 
disposal in building good creation operators is smearing. The best means of understanding the low-energy degrees 
of freedom of confinement must thus come from a smearing method that can address all these features; it must be 
numerically accessible, do a good job of projecting onto the lowest energy states and it must facilitate easy evaluation 
of correlation functions involving a broad set of creation operators, including those exciting multi-hadron states. 
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A new generation of experiments devoted to hadron spectroscopy, including GlueX at Jefferson Lab. PANDA at 
GSI/FA1R and BES 111 intend to make measurements with unprecedented precision and in previously unexplored 
mass ranges and quantum-number sectors. The main aim of the Hadron Spectrum Collaboration is to extract from 
realistic lattice simulations predictions for masses, decay properties and relevant matrix elements in these domains. 

Spectroscopy investigations are almost as old as non-perturbative Monte Carlo calculations on the lattice 0, [|| . 
As the first studies progressed, it was quickly realised that operators creating hadrons which are constructed directly 
from the fields in the lattice lagrangian have significant overlap with a large tower of states and that extracting the 
lightest elements of the spectrum was difficult. The solution was to build operators from "smeared" fields, where some 
linear operator was applied first to the quark degrees of freedom on the appropriate time-slice before the creation 
operator was formed. The aim of constructing a smeared field was to reduce the component of the fields close to the 
cut-off substantially, since these modes do not contribute significantly to long-range correlation functions. Initially, 
the techniques considered smeared fields that were not gauge covariant. Gauge-covariant schemes were introduced 
soon after 0,||. In particular, Ref. 0] introduced the Jacobi smearing algorithm, in which a lattice approximation to 
the three-dimensional gauge-covariant laplacian is constructed and applied iteratively to the field. This acts naturally 
as a long-wavelength filter on the modes constituting the quark field. 

In this paper, a new method for smearing quark fields is described and tested in large-scale realistic simulations. 
In Sec. HH the formulation of the new algorithm is presented, including detail on its application to measurements 
of hadron two-point correlation functions and three-point functions where an arbitrary current is inserted. Sec. IIIII 
presents the results of initial tests of the effectiveness of the method, and Sec. II VI discusses practical issues that arise 
with the method and outlines some future research directions. 



II. OPERATOR CONSTRUCTION 



The energy of an eigenstate of the hamiltonian of a quantum field theory can be determined by computing the 
correlation function between creation and annihilation operators x a t Euclidean times t and t'; 

c(t',t) = ( x (t')xHt))- (1) 

Inserting a complete set of eigenstates of the hamiltonian, such that H\k) = Ek\k), this correlation function decom- 
poses into a sum of contributions from all states in the spectrum with the same quantum numbers as the source 
operators, 

C(t\t) = Y,\(x\k)\ 2 e~ E ^'~ t l (2) 

k 

In order to measure energies of low-lying states, it it crucial to construct operators that overlap predominantly with 
these light modes. This exposes the asymptotic behaviour of the correlator at earlier time-separations and enables more 
statistically accurate determinations of energies. There is a good deal of freedom in the construction of appropriate 
operators. Once the constraints of symmetries and temporal localisation are imposed, any function of the fields in 
the path integral can be used in principle. 

Smearing is a well-established means of defining an initial step in the construction of creation operators. Rather 
than applying a creation operator to the fields directly, a smoothing function is applied first. The function should 
preserve as many symmetries as possible while effectively removing the presence of short-range modes, which make 
an insignificant contribution to the low-energy correlation function. In contemporary simulations, a popular form of 
this operation in theories where fermion fields couple to gauge bosons is the Jacobi smearing method |8(. 

Gauge-covariant quark smcarings based on the lattice Laplacian start with the simplest representation of the 
second-order three-dimensional differential operator 
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-V 2 xy (t) = 6S xy - (Uj(x, t)S x+j , y +U]{x- j, t)*x-j,„) , (3) 

where the gauge fields, U may be constructed from an appropriate covariant gauge- field-smearing algorithm @. After 
defining the Laplace operator, a simple smearing can be written 

, (() _( 1 + £?MY*, (4) 
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where er and n a are tunable parameters that can be used to optimise projection onto the states under investigation. 
For large n a , this approximates the exponential of cV 2 , i.e. 

lira J a , n „{t) =cxp(aV 2 (t)). (5) 

n a — >oc 

The resulting exponential suppression of higher eigenmodes of the lattice Laplace operator means only a small number 
of the lowest modes contribute substantially to J. 

This observation suggests the smearing operator can be approximated by forming an eigenvector representation, 
truncated to the lowest modes. Since the smearing operators used in lattice calculations are approximated well by 
low-rank constructions, a definition of smearing can be chosen to enforce this more absolutely. Let Vm be the vector 
space of scalar fields charged under the fundamental representation of the gauge group on a particular time-slice. 
Vm has rank M = N c x N x x N y x N z where N c is the number of colours, and N x ,N y and N z are the extents of 
the lattice in the three spatial directions. Now define smearing to be a well-chosen operator of rank TV <C M . This 
class of operators will be called "distillation" operators. There is a substantial benefit to doing this: if the rank of 
the operator is sufficiently small, all elements of the propagation matrix from this space can be constructed at an 
affordable computational cost (which is proportional to the rank of the space, TV). Consequently, correlation functions 
involving arbitrarily intricate hadron creation operators can be measured with a fixed inversion overhead. 

Define the distillation operator on time-slice t as a product of an M x TV matrix and its hermitian conjugate: 

N 

□ (t) = V{t)V\t) U xy (t) = X>< fc > (6) 

k=l 

The fc th column of V(t) contains the k th eigenvector of V 2 evaluated on the background of the spatial gauge fields 
of time-slice t, once the eigenvectors have been sorted by eigenvalue. This is the projection operator into Vat, the 
subspacc spanned by these eigenmodes, so D 2 = □. When the number of eigenvectors included is the same as the 
dimension of Vm, i.e. TV = M, the distillation operator becomes the identity, and fields acted upon are unsmeared. 

The Laplace operator inherits many symmetries of the vacuum. It transforms like a scalar under rotations, is 
covariant under gauge transformation and is parity and charge conjugation invariant. If the action of one of these 
symmetries on the Laplace operator maps V 2 onto V 2 , then there is a unitary transformation, R on Vm such that 

RV 2 R*< = V 2 . (7) 

This implies that if v is an eigenvector of V 2 then the eigenvectors of V 2 are Rv. Considering the definition of the 
distillation operator given in Eq. [51 the transformed operator must then obey 

R □ ijt = □, (8) 

and so correlation functions constructed using distilled fields have the same symmetry properties on the lattice as 
those constructed using Laplacian smearing methods. 



A. Meson two-point correlation functions 

Consider the momentum-projected creation and annihilation operators of an isovector meson, uT A d and dT B u, 
where F acts in spin and color as well as coordinate space. Applying the distillation operator □ onto each quark field, 
the creation operator at three- momentum p is written as 

xUf,t) = u x (t)n xy (t) ■ e -*-»r£(t) • □ !B (t)4(t), (9) 

where there is an implied summation over repeated spatial indices. In a shorthand notation the correlation function 
can be written as 

C M \t',t) = (d(t')n(t')T B (t')n(t')u(t') ■ u(t)a(t)T A (t)n(t)d(t)). (10) 

After evaluating the quark field path-integral and inserting the outer-product definition of the distillation operator □ 
from Eq. [HI the correlator can be written 

$ B (t')r(i',t)$ A (t)r(t,t')l ! (11) 



C M \t',t)=Tt 
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where 

= vHt) [r A (t)] a0 v(t) = v\t)v A (t)v(t)s^ (12) 



and 



r a p(t',t) = V\t')M-i(t',t)V(t), (13) 



with M the lattice representation of the Dirac operator and where the quark spin indices, a, (3 of $ and t have been 
explicitly written. $ has a well-defined momentum, while there is no explicit momentum projection in the definition 
of r. Often, $ can be decomposed into terms that act only within coordinate and color space V A and only within 
spin space S A . Note that and r are square matrices of dimension N x N a where N a is the number of components 
in a lattice Dirac spinor. Therefore it requires just N x N a operations of the inverse of the fermion matrix on a vector 
in order to compute all elements of r, the "perambulator" . Notice also that the choice of source and sink operators 
is entirely independent of the computation of t; any source and sink operators can be correlated a posteriori once all 
elements of the r matrix have been computed and stored. The method straightforwardly extends to the determination 
of correlation functions for mesons composed of different, non-degenerate quark flavors. 

In the determination of an isoscalar meson correlation function, evaluation of disconnected terms is required. The 
disconnected diagram can be similarly represented once distilled fields arc used in creation and annihilation operators. 
Such a term would comprise two separate traces over the distillation space: 

C^' disc \t',t) = Trh A (t)T(t,t)}Trh B (t')^(t' ,t')\. (14) 

An exact determination of this expression requires computing r(t,t) for all time-slices t. 

Since distilled single-particle operators can be projected onto definite momentum at both source and sink, multi- 
meson correlators can also be constructed using creation and annihilation operators of the form 

XMM(q=Pi -P%]t) = XM(pi,i)XM(-P2,t), (15) 

where p\ and P2 are the three-momenta of the single particle operators in Eq. [9l After integration over the quark 
fields, the resulting diagrams can again be computed by taking traces over products of construction operators in the 
distillation space with the perambulators. The precise structure of these functions depends on the quark flavor content 
of the source and sink operators, so details are not presented here. As before, these matrices are small compared to 
the dimension of the space of quark fields, and once the quark propagation is encoded in the perambulator, these 
multi-hadron correlation function can be computed. 



B. Baryon two-point correlation functions 



The factorization technique can also be applied to baryons. To illustrate the concepts involved, consider just the 
isospin-1/2 sector although the technique generalizes naturally to other baryon multiplcts. An annihilation operator 
involving displacements as well as coefficients in spin space is written: 

XB (t) = e^S^^iV.Ddr^aut^aur^t), (16) 

where the color indices of the quark fields acted upon by the displacement operators T>i are contracted with the 
antisymmetric tensor, and repeated spin indices are summed. After integration over quark fields, the correlation 
function 

C ( B ) (t',t) = -{ X B(t')XB(t)) (17) 

can be factored into perambulator terms, Eq. [13] as well as creation and annihilation operators, where the latter can 
be written as 

= (Viv (l) ) a (V 2 v^)" (V 3 v^y (t) S aia2a3 . (18) 
The spin terms factor from the antisymmetric contraction of the vectors v. 
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There are two terms in the resulting correlation function that involve tensor contractions of the creation and 
annihilation operators with the perambulators 

Cf fa, Tu, Tu] (t\ t) = 4>(^ fc ) (t')T { ; & (t', t)T^ (t\ t)T^ (i', (t) 

_$( i ^*)(t')r^(t' J t)TW' S )(^ I t)T^(t' > t)*^'^*(t) ) (19) 

where there is an implicit tensor contraction over the internal spin indices, and where the quark labels d and u are 
used to denote the corresponding flavors of the quark perambulator terms in Eq. 1131 

Similar to the meson case, the choice of source and sink operators is independent of the computation of the 
perambulators and t u . The computation of these matrices can also be shared with computation of meson correlators. 
A baryon correlation function can be evaluated a posteriori using the contractions of the vectors with displacements 
in Eq. (TS] which can also be shared among the source and sink operators. These contractions do not involve spin 
components, thus making the storage of <I> manageable. 

C. Meson three-point correlation functions 

A generic meson three-point function is written 

C {3) (t f ,t,ti) = ((dnr B Du)(t f ) ■ (uTu)(t) ■ (uUT A nd)(U)Y (20) 

where there is no smearing in the operator on timeslice t. The completely connected Wick contraction can be expressed 
as 

(21) 

where the "generalized perambulator" is defined by 

J a p(tf,t,U) = V\t f ) {M-\t f ,t)Y{t)M-\t,U)) ap V{My (22) 

In general, a disconnected term may also appear; this can be factorised as the product of a two-point function 
(evaluated according to the algorithm of Sec. Ill A|) and a disconnected insertion. Since the disconnected trace does 
not involve distilled fields, the evaluation of this insertion is not discussed further at this point. To compute the 
connected three-point correlation function, there arc two distinct sets of inversions needed. The first is from the 
distillation vectors at the source time-slice t, and the second from the sink time-slice t f . The generalised perambulator 
of Eq. [5Uis then constructed by contracting the two solutions from these two sets of inversions. The current insertion 
is encoded in the choice of operator T. Any momentum insertion at the current operator involves a Fourier transform 
on each time-slice t. 



D. Baryon three-point correlation functions 

Baryon three-point correlation functions can also be expressed in terms of the generalized perambulators of Eq. 1221 
Consider a three-point correlation function where a baryon is created on time-slice U, then is acted on by a current 
operator on time-slice t and is subsequently annihilated at tf : 

C<P(t f ,t,U) = -(B(t f )- (uTu)(t) -B{t t )). (23) 

After quark-field integration, this can be written as a sum of both a connected three-point contribution and the 
product of a disconnected insertion and a two-point correlator. As before, the discussion here is restricted to isospin- 
1/2 light baryons although the technique is quite general. Note that the disconnected insertion involves unsmeared 
fields and must be evaluated through some other technique. The two-point contribution follows from Sec. Ill Bl For 
an up-quark insertion in the correlation function 

C(3) (x x -i \ _ abc abeq q 
conn, v 1- / ' l i L i) fc fc °aia20 3 '-'01x0120:3 

((CDiadr ai (V 2 au) b a2 {V 3 n u Y a3 )(t f ) (24) 
• (uTu)(t) ■ {(da^l^uav^mv^iu)), 



C^U(tf,t,U) =Tr 9 B {t f ) J(t/ ) i,ti)* A (ti)7fer t (t/ ) ti)7fe 
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there are four Wick contractions after quark-field integration. The up-quark generalized perambulator from Eq. [22] is 
denoted as 

Mtf,t,ti) = 7t(t / )M- 1 (t / ,t)r(t)M- 1 (t,^)F(t i ) (25) 

and similarly for the down quark. Since the quarks in the insertion come in pairs, the connected three-point correlator 
can be written as the sum of two two-point correlation functions with each up-quark perambulator substituted with 
the corresponding generalized perambulator in turn; 

CV>{t f ,t, U) = [t d , Tu , Tu){t) + C< a > [td,tv, ru](t). (26) 
The down-quark contribution can be written as 

Cf ] (tf,t,U) = C^[T D ,Tu,Tu]{t). (27) 

This kind of construction is a generic feature of the three-point functions. 

As in the previous two- and three-point correlator examples, the choice of source and sink operators is independent 
of the perambulator and generalized perambulator computations. These latter generalized perambulators can be 
shared among meson three-point correlator calculations. In addition, for degenerate quarks, only one generalized 
perambulator need be computed, independent of the flavor of the quark line where the current insertion is to occur. 
Thus, the overall computation of distinct three-point correlation functions has a manageable computational cost. 

III. TEST RESULTS 

To test the efficacy of the method, a first study of meson and baryon two-point correlators was performed. In 
addition, the method is demonstrated for a two-meson system. The background gauge fields were generated by the 
Hadron Spectrum Collaboration [HI EH- They are 16 3 x 128 and 20 3 x 128 anisotropic lattices, with the ratio of 
spatial to temporal lattice spacings a s /a t = 3.5 and 2+1 flavors of dynamical quarks. The bare light-quark mass is 
atm q = —0.0840, and when the f2 baryon mass is used to set the scale, the corresponding pion mass is 383 MeV. The 
valence-quark action is the same as the light-dynamical-quark action. 

The numerical techniques for the iterative solution of the sparse linear system for the Sheikholeslami-Wohlert 
improved discretisation of the Dirac equation are described in Ref. In particular, the EigCG inverter [12] was 

used for the solution of these linear systems. This method uses eigenvector deflation to accelerate convergence. The 
first few inversion calls are needed to construct a sufficient basis, and subsequent inversions are accelerated by a factor 
of about three or more compared to the undeflated solver. For the lattices considered here, this asymptotic rate is 
found for N > 4. 



A. The profile of the distillation operator 

The distillation operator is written as a low-rank projection operator into a vector space of smooth fields. As such, 
there is little physical intuition to suggest it resembles more traditional smearing algorithms in generating a field that 
is roughly localized in some confining region while still remaining smooth. To examine the properties of the operator, 
its spatial distribution was computed. Defining 

m (28) 

means "J is a gauge-invariant measure of the degree to which a field is smeared by the application of the distillation 
operator. This observable was measured for various values of N and for a range of different off-set vectors r on a 
small ensemble of the 20 3 spatial lattices and the results are presented in Figure [TJ The data clearly shows that for 
N between 8 and 64 the distillation operator is very well described by a gaussian wavefunction. As N increases, 
so the radius of this distribution reduces. This behaviour is to be anticipated since when N = M the distillation 
operator is the identity and has no spatial extent at all. The measurements presented in Fig. [T]also clearly illustrate 
the distillation operator has good rotational symmetry, indicating it is only very weakly influenced by dynamics at 
the lattice cut-off. The artifacts seen at larger values of r in the very broad TV = 8 distribution arise solely from the 
finite lattice volume used in these measurements. 
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FIG. 1: The spatial distribution of the distillation operator constructed from the lowest eigenvalues of the Laplace operator 
on the lattice. Data is from a 20 3 spatial volume and statistical fluctuations are much smaller than plot symbols. The dashed 
lines are fits to a gaussian distribution. 
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TABLE I: Derivative-based meson operators used in these tests. V is the one-link discretised covariant derivative 



B. Meson correlation functions 



In this section, connected meson correlation functions formed using the derivative-based operators presented in 
Refs. [l3|, EH are shown, with particular focus on the T x channel. Fig [2] shows plots of the effective mass, defined 

by 

m ^ t)= .uJ^m (29) 



st \ C{t) 

with St = bat, for four diagonal correlators while varying the number of distillation vectors. Also shown are the corre- 
sponding effective masses computed using the traditional point-to-all propagator construction with both unsmeared 
and Jacobi-smeared sources. The four operators used are defined in TableU In all cases, gauge links in the operators 
are stout-smeared @. 90 configurations of the 16 3 x 128 lattice were used for this comparison. 

A variational analysis [HI, [Ty| of matrices of two-point correlation functions has proven to be a powerful method 
in extracting the spectrum and matrix elements of excited states [l7| . Consequently, cost estimates in terms of the 
computing time required to form such a matrix of correlators and not just any single correlator should be considered. 
As a measure of the relative computational demands of distillation versus the traditional point-to-all approach, note 
that in the point-to-all approach, ten forward propagators, generated by applying I, ^%= x ,y,z, ViVjj^j to a point 
source at the origin are required to produce the set of T x correlation functions defined above. Each of these sources 
then requires twelve inversions of the Dirac matrix (N c x N a ). To compare with the distillation methods, note that 
approximately the same total inversion cost arises from using 32 distillation vectors since 

(Ar srcs = 10) x (A^ c = 3) x (A^ = 4) = 120 w (N = 32) x (A^ = 4) = 128. (30) 

In the correlation function for the ip"fiip operator, it is clearly seen that increasing the number of vectors causes 
the effective mass to approach that of the unsmeared point-to-all correlation function, since ^fe=i v^v^^ = 1. The 



8 





FIG. 2: Effective masses (five timesiice shift) for diagonal correlators in the T x channel computed on 16 3 x 128 lattices. The 
operator construction is shown in table [I] Correlators are constructed using distilled quark fields with N 6 {16, 24, 32, 48, 64} 
vectors and using point-to-all methods with standard gauge-invariant quark smearing (SS, brown) and without (PP, black). 



statistical variance is observed to decrease rapidly as the number of vectors is increased. For this range of N it is 
always smaller than the point-to-all correlators. A similar but less rapid convergence of the signal is observed in the 
«o x Vtj channel. 

On the other hand, the it x signal appears to be essentially unchanged over this range of vector number but 
the noise reduces considerably with increasing vector number. The p x ESj^ channel shows similar behaviour. Note 
that these operators sample non-trivial spatial dependence, and this may require a larger space of eigenvectors. The 
rapid scaling of noise with number of vectors can be quantified as follows: consider the correlator noise-to-signal ratio 
r on a timesiice as a function of the number of vectors used, N and model this behaviour with 

r = a+—. (31) 

Fitting the parameters to best reproduce the Monte Carlo data from 450 configurations and averaging over many 
time-slices gives the following estimates for the scaling exponent 

p( 7< ) = 0.8(1); p(o x V Tl ) = 1.1(2); p(tt x B Tl ) = 1.3(3); p(p x lD> Tl ) = 1.8(2). 

In all cases, the noise reduction is better than the behaviour expected from simple statistical scaling, p = 0.5. 

As a realistic test of the efficacy of distillation in the extraction of an excited-state spectrum, a matrix of correlation 
functions using twelve operators with quantum numbers was constructed. In Fig. [3] the low-lying spectrum 

extracted as a function of number of vectors used is shown. Clearly, even with a relatively small number of vectors 
this method can be used to construct two-point correlation matrices from which the spectrum of low-lying states can 
be extracted with some precision. The masses, up through the fourth excited state, become more consistent with 
increasing number of vectors. 

Now consider the issue of the scaling of this method with increasing lattice volume. Fig. 2] shows the effective 
masses of correlators computed on both 16 3 and 20 3 lattices, the spatial-volume ratio here being very close to 2 
(|p- = 1.95). It is evident that the correlation function with N vectors on the 16 3 volume is essentially the same as 
the corresponding case with 2N vectors on the 20 3 lattice. Note that although the number of inversions for the same 
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FIG. 3: T 1 variational mass spectrum as a function of number of eigenvectors used. The error bars are purely statistical and 
do not take account of possible systematic errors due to fitting. 



signal is doubled when the lattice volume increases from 16 3 to 20 3 , there is a reduction in the noise of approximately 
a factor of two. 

A strikingly simple explanation for this effect comes through consideration of the eigenvalues belonging to these 
eigenvectors of the Laplacian operator, S/ 2 v^ = —XkV^ k '. Fig. [5ja) shows the eigenvalues measured on 16 3 and 20 3 
lattices. One could consider choosing the rank of the distillation operator by using all eigenvectors with eigenvalue 
below a certain cut-off. It is clear from the plot that on a larger volume this will require more eigenvectors. This is 
displayed more directly in Fig. [HJb) where the number of eigenvectors with eigenvalue less than a cut-off is plotted. 
Approximately twice as many eigenvectors in the 20 3 case are required to have the same cut-off as in the 16 3 case. 
It appears from this analysis that the meson correlator signal is effectively determined by the cut-off in eigenvalue 
used to construct the distillation space. This is equivalent to using a Heaviside approximation to smearing: J\(t) = 

e(A + v 2 (t)) = Efci e(A - A fc )« (fc) (t)« (fc)t (t) = Ekli « (fc) (*y fe)t (*)- 



C. Baryon correlation functions 

Baryon correlation functions were evaluated using the displaced-quark operators described in Refs. (l8l . [l9| and 
employed in spectrum studies Refs. [20L l2lj. The distillation method is demonstrated on four nucleon G\ g operators, 
two of which are local to a single site with the remaining two having a singly displaced quark field. The corresponding 
4x4 matrix of correlators was computed on 316 configurations of the 16 3 x 128 lattice. Fig.[5]shows the effective-mass 
plots for the diagonal correlators. A clear trend is observed, similar to that seen for the mesons, of the correlator 
having larger excited-state contribution but smaller statistical fluctuations as the number of eigenvectors is increased. 
Modelling the correlator noise-to-signal ratio with Eq. |3"T1 gives best-fit exponents p ~ 1-1(2), showing again that 
increasing the number of vectors decreases the noise considerably faster than simple statistical scaling. The variational 
method [HI, [l6| was used to extract the masses of the lowest two states in the G\ g spectrum. Fig. [7] shows the 
dependence of the effective masses of these states on the number of eigenvectors used to form the distillation operator. 
Consistent masses are found, with an increase in statistical precision as the number of vectors is increased. Fig. [H] 
shows the diagonal effective masses computed on 16 3 and 20 3 lattices. The similarity between the signals found before 
and after doubling the number of vectors is not as obvious here as in the meson case. 
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Operator construction shown in table [I] 
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D. Multi-meson correlation functions 



One of the main motivations for defining distillation methods is to enable efficient evaluation of correlation functions 
of creation operators which resemble multi-hadron states where the component hadrons have definite momentum. 
Reliable and efficient means of including these states in the operator basis are important for addressing scattering 
and resonance properties. The simplest case was tested, which comprises the isospin-two two-pion state formed using 
Eq. ll5l with the single particle operators M(p Te \) constructed using the local operator u^^d. For this example, consider 
the two-pion operator 



XMM(\Piel\ 2 ,t) = Xu{k,t)XM{-k,t) 



(32) 
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with total momentum zero, and the single-particle operators to have relative momentum p re \. The sum is over the set of 
momentum vectors lZ(p ro \) related to p TC i by lattice rotations, to create a state that transforms trivially under rotations 
and so the operator forms an Ai irreducible representation of the cubic group. In this test, all values of |p r oi| between 
and 2x^ were considered, corresponding to p rcl = ^-n with ft € {(0,0,0), (1,0,0), (1, 1,0), (1, 1, 1), (2,0,0)}. 

A 5 x 5 correlation matrix is formed between all of these two-pion operators which is diagonalized by solving the 
generalized eigenvalue problem. The effective mass of the lowest three principal correlation functions measured on 
100 configuration is shown in Fig. [9l Clear indication of a signal for the second excited state is seen, even with 
this relatively low level of statistics. Note that when using the conventional point-to-all method, the source operator 
cannot be projected onto definite momentum so the correlation matrix described above cannot be constructed in 
practice. 
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IV. DISCUSSION AND FUTURE DIRECTIONS 



The results of Sec. [TTT] demonstrate clearly that the light hadrons can be studied using operators that excite only a 
very small number of extended degrees of freedom on each time-slice. The restricted set of modes needed to create 
a hadron means all elements of quark propagation needed to construct the correlation function can be measured 
directly. This does not preclude using stochastic methods [|| to reduce computational costs. 

In the straightforward recipe described in this paper, the distillation space was constructed from the lowest eigen- 
vectors of the lattice three-dimensional Laplace operator on a time-slice. This definition is not unique and more 
experimentation with different choices of distillation space might prove useful. Also, the operator used in this work is 
diagonal in the quark spin index. Alternative choices of the distillation space might include non-trivial spin structure, 
using for example the vector space spanned by the lowest modes of the three-dimensional Dirac operator. 
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Recall that in the simplest implementation tested here when N = M, the distillation operator becomes the identity 
and no smearing is performed. This is clearly an undesirable feature that is circumvented by the obvious variation on 
the method to include different weights for the eigenmodes. This more general distillation operator is then 

N 

J(t)=J2v {k) (t)f(\ k (t))v^(t). (33) 

fc=i 

Including weight functions on the eigenvectors can be used to suppress more rapidly fluctuating contributions when 
many eigenvectors are employed. Using exponential weights leads to increasingly accurate approximations to the 
conventional smearing algorithm of Eq. 0J 

No increase in the number of eigenvectors would be expected as the continuum limit is taken at a fixed physical 
volume since, as is seen in the data of Sec. IIII A[ the smooth distilled modes are almost completely decoupled from 
the cut-off dynamics. This expectation is yet to be tested numerically. One substantial possible drawback with the 
method is the need to increase the rank of the distillation operator as the three-dimensional lattice volume increases. 
The corresponding linear increase in the cost of solving the Dirac equation leads to an algorithm for spectroscopy that 
scales like 0{V 2 ). It appears from the results we have presented that with the required increase in rank with increasing 
volume one also obtains a decrease in statistical fluctuations. The possible solution to the issue of volume scaling may 
be to combine the method with a suitably designed stochastic estimation scheme. This is under investigation. 

Note that the method does not give direct access to all elements of the quark propagator, only those relevant to 
low-energy spectroscopy. This would be a limitation in calculations involving isoscalar operator insertions, where 
disconnected loops with all momentum components of the quark field are needed. A well-known example of such a 
computation would be the evaluation of the strangeness content of the nucleon. 

V. SUMMARY 

In this paper, a simple new quark-field smearing has been proposed and tested in numerical investigations using 
QCD gauge fields with N j = 2 + 1 dynamical flavors. There is substantial freedom in the construction of creation 
operators for hadrons. This has been exploited in this work to define an algorithm that applies a low-rank projection 
operator to the quark fields of the path integral on each time-slice in order to define smooth modes for subsequent 
operator construction. This makes affordable the evaluation of all information about quark propagation required 
to measure distilled hadronic correlation functions. Once the propagation matrix has been evaluated, correlations 
between arbitrary choices of creation and annihilation operators can be determined without requiring any further 
inversions. The theoretical framework needed to extend the scope of measurements using this technology to include 
matrix elements was described briefly. 

Some initial tests of the method demonstrate that restricting quark fields to lie in a very low-dimensional space 
of smooth fields does not substantially degrade the quality of Monte Carlo evaluations of correlation functions and 
in many cases, statistical accuracy is enhanced. The most substantial drawback, still to be overcome, is the growth 
in cost of the algorithm as the lattice volume in physical units increases. Since the method is concerned with long- 
distance physics, no extra difficulties in approaching the continuum limit are anticipated, but this would need to be 
confirmed by practical testing. 

The method is currently in its infancy and a number of research directions are outlined in this paper. The primary 
focus of the development efforts of this collaboration is to include stochastic evaluation of correlation functions in an 
optimal way and results will appear soon. The collaboration will use the techniques outlined in this paper in a range 
of spectroscopy determinations in the near future. 
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